# This is a python file created by AJB in Sept 2025 to create a plot for the paper
#
#AJB here is the command: python3.11 ./ajb_plot-mars-visc.py
#here is the colour page:
#https://matplotlib.org/2.0.2/examples/color/colormaps_reference.html
import matplotlib.pyplot as plt
import numpy as np

import matplotlib as mpl
from matplotlib.colors import LinearSegmentedColormap
#fig, ax = plt.subplots()             # Create a figure containing a single Axes.
#ax.plot([1, 2, 3, 4], [1, 4, 2, 3])  # Plot some data on the Axes.
#plt.show()                           # Show the figure.

np.random.seed(19680801)

z = [219.14,
284.72,
339.98,
380.69,
152.97,
141.6,
168.71,
200.98,
239.65,
283.63,
179.54,
211.35,
250.29,
299.66,
367.6,
347.57,
186.9,
231.11,
254.27,
156.23,
185.4,
221.31,
267.11,
163.7,
156.9,
199.26,
151.48,
145.88,
176.78,
214.37,
318.51,
314.46,
283.89,
260.55,
259.1]

#B's are y's
y = [2571.6,
2571.6,
2571.6,
2571.6,
6171.6,
6671.6,
6671.6,
6671.6,
6671.6,
6671.6,
6171.6,
6171.6,
6171.6,
6171.6,
3571.6,
4571.6,
4571.6,
4571.6,
3571.6,
5161.6,
5161.6,
5171.6,
5171.6,
2571.6,
3571.6,
3571.6,
4571.6,
5571.6,
5571.6,
5571.6,
5571.6,
3671.6,
4671.6,
5571.6,
5671.6]

#C's are x's
x = [606.2,
506.2,
406.2,
306.2,
706.2,
706.2,
606.2,
506.2,
406.2,
306.2,
606.2,
506.2,
406.2,
306.2,
306.2,
306.2,
606.2,
506.2,
506.2,
706.2,
606.2,
506.2,
406.2,
706.2,
706.2,
606.2,
706.2,
706.2,
606.2,
506.2,
306.2,
406.2,
406.2,
406.2,
406.2]


#numpy_array = np.array([x, y, z])
fig, ax = plt.subplots()
#for color in ['tab:blue', 'tab:orange', 'tab:green']:
for i in range(2,len(z)):
    n = len(z)
    maxval = max(z)
#    x, y = np.random.rand(2, n)
#    scale = 200.0 * np.random.rand(n)
    scale = 10.0 * z[i]
    ax.scatter(x[i], y[i], c=z[i]*10.0/maxval, s=scale, label=i,
               alpha=0.3, edgecolors='none', cmap='RdYlBu')

    
#ax.legend()

ax.set_title("AJB paper Figure")
ax.grid(True)

plt.show()


# let's try again:
#next add x, y axis names and colorbar name and change colour style
#plt.scatter(x, y, c=z, cmap='viridis')
plt.scatter(x, y, c=z, cmap='rainbow', s=4000, marker='s')
plt.xlabel("C (K)", size=12)
plt.ylabel("B (kJ/mol)", size=12)

plt.grid(True)
cbar = plt.colorbar()
cbar.set_label("Flow distance in km", labelpad=15, rotation=270)
#let's try to add the Seitah pts on top
plt.scatter(406.2, 5671.60, c=259.1, cmap='gray', s=25, marker="P")
plt.scatter(610.96, 5442.27, c=181.86, cmap='gray', s=25, marker="P")
plt.scatter(618.98, 5776.58, c=178.06, cmap='gray', s=25, marker="P")
plt.scatter(735.00, 3164.38, c=145.73, cmap='gray', s=25, marker="P")
plt.scatter(667.36, 4200.14, c=170.06, cmap='gray', s=25, marker="P")
plt.scatter(714.82, 3292.56, c=158.56, cmap='gray', s=25, marker="P")
plt.scatter(588.28, 4888.46, c=192.22, cmap='gray', s=25, marker="P")
plt.scatter(600.81, 5295.16, c=189.16, cmap='gray', s=25, marker="P")
plt.scatter(601.31, 5279.28, c=186.16, cmap='gray', s=25, marker="P")
plt.scatter(615.45, 5946.85, c=177.38, cmap='gray', s=25, marker="P")
#added in a label next change the size, 
#put in other points, change to big cross
plt.annotate("Séítah", (x[len(x)-1]+5,y[len(y)-1]+5), size=8, color='indianred')
plt.annotate("Artuby", (610.96+5, 5442.27+5), size=8, color='salmon')
plt.annotate("Máaz", (618.98+5, 5776.58+5), size=8, color='tomato')
plt.annotate("BAMb", (735.00+5, 3164.38+5), size=8, color='dodgerblue')
plt.annotate("IAMb", (667.36+5, 4200.14+5), size=8, color='royalblue')
plt.annotate("AdMb", (714.82+5, 3292.56+5), size=8, color='steelblue')
plt.annotate("PIXL Dourbes", (588.28+5, 4888.46+5), size=8, color='olivedrab')
plt.annotate("65-3", (600.81+5, 5295.16+5), size=8, color='forestgreen')
plt.annotate("65-5", (601.31+5, 5279.28+5), size=8, color='olive')
plt.annotate("Di9010", (615.45+5, 5946.85+5), size=8, color='darkslategray')
plt.title("Flow distance for variable viscosity lavas", size=16)
plt.show()


# let's try to make the colours map work here
x = np.arange(0, np.pi, 0.1)
y = np.arange(0, 2 * np.pi, 0.1)
X, Y = np.meshgrid(x, y)
Z = np.cos(X) * np.sin(Y) * 10
colors = [(1, 0, 0), (0, 1, 0), (0, 0, 1)]  # R -> G -> B
n_bins = [3, 6, 10, 100]  # Discretizes the interpolation into bins
cmap_name = 'my_list'
fig, axs = plt.subplots(2, 2, figsize=(6, 9))
fig.subplots_adjust(left=0.02, bottom=0.06, right=0.95, top=0.94, wspace=0.05)
for n_bin, ax in zip(n_bins, axs.flat):
    # Create the colormap
    cmap = LinearSegmentedColormap.from_list(cmap_name, colors, N=n_bin)
    # Fewer bins will result in "coarser" colomap interpolation
    im = ax.imshow(Z, origin='lower', cmap=cmap)
    ax.set_title("N bins: %s" % n_bin)
    fig.colorbar(im, ax=ax)

#ax.legend()
ax.grid(True)

plt.show()



#now let's try a line plot for Figure 3b:, this is the host (y) and parasite (X) population plot
#from mpl_toolkits.axes_grid1 import host_subplot
from matplotlib.collections import EventCollection
from matplotlib.patches import Rectangle
#ajb this creates a twin x vs y axis
#host = host_subplot(111)
#par = host.twinx()
fig, ax = plt.subplots()

ax.set_xlabel(r'relative viscosity at 1400K $\frac{viscosity \nu}{{\nu _{Di_{9O}An_{10}}}}$')
ax.set_ylabel(r"Depth of crystals relative to lake depth - $\frac{h_f}{H}$")
ax.add_patch(Rectangle((20, 0), 80, 0.02, color='grey', alpha=0.5))
#par.set_ylabel("Temperature")
#this will return the Line2D objects p1, p2:
p1, = ax.plot([1, 10, 100, 10000, 1000000],[0.52, 0.45, 0.41, 0.4, 0.395], label="Bottom cooled", color="blue", marker='o', linestyle='-')
p2, = ax.plot([1, 10, 100, 10000, 1000000],[0.32, 0.18, 0.1, 0.02, 0.01], label="Bottom insulated", color="red", marker='*', linestyle='--')

#set the y axis limits based on Worster
plt.ylim(0, 0.6)
plt.xlim(0.01, 100000)

ax.legend(labelcolor="linecolor")
ax.legend(loc='center right')

#ax.yaxis.label.set_color(p1.get_color())
#ax.yaxis.label.set_color(p2.get_color())

#line, = ax.plot(xitem, yitem)

x_data = p1.get_xdata()
y_data = p1.get_ydata()

p1.set_color(p1.get_color())
p2.set_color(p2.get_color())

visc_data = [0.06,0.08,0.38,0.392,0.49,0.68,0.69,0.77,0.94,1]

#lets put in the labels (annotations) for the viscosities:
# Create an EventCollection for x-axis events
xevents = EventCollection(visc_data, color='red', linelength=0.02, lineoffset=0)
# Label a specific event, for example, the 3rd event
event_index_to_label = 2
line_opacity = 0.2
line_shrink = 0.005
#ax.annotate(f'Séítah {event_index_to_label}', color='red', rotation=45,
ax.annotate(f'BAMb', color='dodgerblue', rotation=70,
            #xy=(x_data[event_index_to_label], 0), # Position at the event on the x-axis
            xy=(0.06, 0), # Position at the event on the x-axis
            xytext=(3, 0.5), # Offset the text
            arrowprops=dict(facecolor='dodgerblue', shrink=line_shrink, alpha=line_opacity))
ax.annotate(f'AdMb', color='steelblue', rotation=70,
            #xy=(x_data[event_index_to_label], 0), # Position at the event on the x-axis
            xy=(0.08, 0), # Position at the event on the x-axis
            xytext=(5, 0.5), # Offset the text
            arrowprops=dict(facecolor='steelblue', shrink=line_shrink, alpha=line_opacity))   
ax.annotate(f'Séítah', color='indianred', rotation=70,
            #xy=(x_data[event_index_to_label], 0), # Position at the event on the x-axis
            xy=(0.38, 0), # Position at the event on the x-axis
            xytext=(8, 0.5), # Offset the text
            arrowprops=dict(facecolor='indianred', shrink=line_shrink, alpha=line_opacity))      
ax.annotate(f'IAMb', color='royalblue', rotation=70,
            #xy=(x_data[event_index_to_label], 0), # Position at the event on the x-axis
            xy=(0.39, 0), # Position at the event on the x-axis
            xytext=(13, 0.5), # Offset the text
            arrowprops=dict(facecolor='royalblue', shrink=line_shrink, alpha=line_opacity)) 
ax.annotate(f'65-3', color='forestgreen', rotation=70,
            #xy=(x_data[event_index_to_label], 0), # Position at the event on the x-axis
            xy=(0.49, 0), # Position at the event on the x-axis
            xytext=(20, 0.5), # Offset the text
            arrowprops=dict(facecolor='forestgreen', shrink=line_shrink, alpha=line_opacity))      
ax.annotate(f'Dourbes', color='olivedrab', rotation=70,
            #xy=(x_data[event_index_to_label], 0), # Position at the event on the x-axis
            xy=(0.68, 0), # Position at the event on the x-axis
            xytext=(30, 0.5), # Offset the text
            arrowprops=dict(facecolor='olivedrab', shrink=line_shrink, alpha=line_opacity))      
ax.annotate(f'65-5', color='olive', rotation=70,
            #xy=(x_data[event_index_to_label], 0), # Position at the event on the x-axis
            xy=(0.69, 0), # Position at the event on the x-axis
            xytext=(50, 0.5), # Offset the text
            arrowprops=dict(facecolor='olive', shrink=line_shrink, alpha=line_opacity))          
ax.annotate(f'Artuby', color='salmon', rotation=70,
            #xy=(x_data[event_index_to_label], 0), # Position at the event on the x-axis
            xy=(0.77, 0), # Position at the event on the x-axis
            xytext=(75, 0.5), # Offset the text
            arrowprops=dict(facecolor='salmon', shrink=line_shrink, alpha=line_opacity))   
ax.annotate(f'Máaz', color='tomato', rotation=70,
            #xy=(x_data[event_index_to_label], 0), # Position at the event on the x-axis
            xy=(0.94, 0), # Position at the event on the x-axis
            xytext=(120, 0.5), # Offset the text
            arrowprops=dict(facecolor='tomato', shrink=line_shrink, alpha=line_opacity))   
ax.annotate(r'$Di_{9O}An_{10}$', color='darkslategray', rotation=70,
            #xy=(x_data[event_index_to_label], 0), # Position at the event on the x-axis
            xy=(1, 0), # Position at the event on the x-axis
            #xytext=(x_data[event_index_to_label] +0, 0.5), # Offset the text
            xytext=(200, 0.5), # Offset the text starts at this x,y
            arrowprops=dict(facecolor='darkslategray', shrink=line_shrink, alpha=line_opacity))
                        
ax.add_collection(xevents)


#print some debug messages:
#print(p1.get_color()) 
#print(x_data)
# Get all Line2D objects from the Axes
all_lines = ax.get_lines()
# You can now access individual lines and their properties
print(f"Number of lines on the Axes: {len(all_lines)}")

print(f"X-data of the first line: {x_data[:5]}...")

for line in all_lines:
    print(f"Line label: {line.get_label()}")



ax.semilogx(x_data, y_data)
# here's another way of doing log scale:
#t = np.arange(0.1, 100000.0, 1)
#ax.xaxis._set_scale('log', base=10)
 
#ax.set_title('Figure 3b')

plt.show()

from matplotlib.collections import EventCollection
# Fixing random state for reproducibility
np.random.seed(19680801)

# create random data
xdata = np.random.random([2, 10])

# split the data into two parts
xdata1 = xdata[0, :]
xdata2 = xdata[1, :]

# sort the data so it makes clean curves
xdata1.sort()
xdata2.sort()

# create some y data points
ydata1 = xdata1 ** 2
ydata2 = 1 - xdata2 ** 3

# plot the data
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.plot(xdata1, ydata1, color='tab:blue')
ax.plot(xdata2, ydata2, color='tab:orange')

# create the events marking the x data points
xevents1 = EventCollection(xdata1, color='tab:blue', linelength=0.05)
xevents2 = EventCollection(xdata2, color='tab:orange', linelength=0.05)

# create the events marking the y data points
yevents1 = EventCollection(ydata1, color='tab:blue', linelength=0.05,
                           orientation='vertical')
yevents2 = EventCollection(ydata2, color='tab:orange', linelength=0.05,
                           orientation='vertical')

# add the events to the axis
ax.add_collection(xevents1)
ax.add_collection(xevents2)
ax.add_collection(yevents1)
ax.add_collection(yevents2)

# set the limits
ax.set_xlim([0, 1])
ax.set_ylim([0, 1])

ax.set_title('Figure 3b')

# display the plot
plt.show()